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A combination of first-principles density functional theory calculations and a search over structures 
predicts the stability of a proton-transfer modification of ammonia monohydrate with space group 
Pi/nmrn. The phase diagram is calculated with the PBE density functional, and the effects of 
a semi-empirical dispersion correction, zero point motion, and finite temperature are investigated. 
Comparison with MP2 and coupled cluster calculations shows that the PBE functional over-stabilizes 
proton transfer phases because too much electronic charge moves with the proton. This over-binding 
is partially corrected by using the PBEO hybrid exchange-correlation functional, which increases the 
enthalpy of P4/nmm by about 0.6 eV per formula unit relative to phase I of ammonia monohydrate 
(AMH-I) and shifts the transition to the proton transfer phase from the PBE pressure of 2.8 GPa to 
about 10 GPa. This is consistent with experiment as proton transfer phases have not been observed 
at pressures up to ~9 GPa, while higher pressures have not yet been explored experimentally. 

PACS numbers: 81.40.Vw,71.15.Nc,31.15.A-,61.66.Fn 



INTRODUCTION 



Ammonia monohydrate (AMH, NH3H2O) exists as at 
least six different crystalline polymorphs over the exper- 
imentally studied range of pressures and temperatures of 
< p < 9 GPa and 170 < T < 295 K± The crystal 
structures of three of these polymorphs have been deter- 
mined: the low pressure phase AMH-I^^ the high pres- 
sure disordered body-centred-cubic phase AMH-VI^ and, 
most recently, a combination of ab initio random struc- 
ture searching (AIRSS)^ and neutron powder diffrac- 
tion data led to the solution of the crystal structure of 
AMH-II^ The structures of AMH-III, IV and V remain 
to be determined. The crystal structures and proper- 
ties of AMH polymorphs (and of the related compound 
ammonia dihydrate, ADH) are of interest to planetary 
scientists due to the likely presence of substantial frac- 
tions of ammonia in ice accreted into the satellites of 
the Gas Giant planets^ilil Whilst the water-rich com- 
pound (ADH) may have a greater abundance than AMH 
at low pressures, it is known that a high-pressure form 
of ADH becomes unstable with respect to a mixture of 
high- pressure AMH and water ice at ^3.5 GPaJ^ Such 
a pressure is relevant to the core of the large icy satellite 
Titan if it undifferentiated (of uniform composition), 12 
as well as during the period of its accretion.— Such a 
pressure may also occur in the icy mantles of fully differ- 



entiated giant icy exoplanets or exomoons. 14 NH3, H2O 
and CH4 are likely to comprise a substantial fraction of 
the interiors of Uranus and Neptune at pressures up to 
600 GPa and temperatures up to 7000 K; 1 — The proper- 
ties of this high p-T molecular mixture are thought to be 
important in the generation of unusual magnetic fields in 
these bodiesJ^ 

Neutron single-crystal diffraction, and the indexing 
and solution of a structure from powder diffraction data 
are non trivial at high pressure. This, combined with 
the tendency towards the formation and/or persistence of 
metastable phases in low-temperature condensed molec- 
ular systems on laboratory timescales, makes it clear that 
there is a central role for the computational prediction of 
equilibrium crystal structures and the determination of 
their physical properties. 

In this work we combine first-principles density func- 
tional theory (DFT) calculations with a random search 
strategy in order to identify new candidate high-pressure 
AMH structures and to compute their stability with re- 
spect to one another and to other known crystalline poly- 
morphs of AMH. Earlier DFT studies of ammonia mono- 
hydrate, ammonia hemihydrate (AHH) and solid ammo- 
nia, have revealed a propensity towards proton trans- 
fer (i.e., formation of an ionic solid) at high pressures. 
Calculations have suggested that AMH-I transforms to 
ammonium hydroxide at ~5 GPa, 17 AHH transforms 
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to ammonium hydroxide ammoniate at ~12 GPa,— and 
solid ammonia transforms to ammonium amide at ^90 
GPa. 19 The ionic solids derived from the two hydrates are 
isosymmetric with their molecular precursors, but there 
is no reason to suppose that other ionic structures might 
not be energetically stable. 

The presence of weak hydrogen bonds and the occur- 
rence of both homo- and hetero-nuclear hydrogen bonds 
in the ammonia hydrates provide a challenge for elec- 
tronic structure methods, particularly with respect to the 
accuracy of exchange-correlation functionals. DFT cal- 
culations using standard functionals such as the Perdew- 
Burke-Ernzcrhof (PBE) generalized gradient approxima- 
tion (GGA) predict hydrogen-bonded molecular phases 
of AMH at low pressures, in agreement with experi- 
ment. We show here, however, that this approach pre- 
dicts molecular AMH phases to be unstable to the forma- 
tion of ionic ammonium hydroxide (NH^:OH~) proton- 
transfer phases at pressures of about 2.8 GPa, although 
no experimental evidence for such phases has been found 
to date, even at pressures up to about 9 GPa. The weak 
van der Waals forces, which are not described by stan- 
dard density functionals such as the PBE-GGA, turn out 
to be important in determining the volumes and relative 
enthalpies of the phases in this system. The zero-point 
(ZP) motion of the H atoms is also important for an 
accurate account of the energetics. We have, moreover, 
found that the most serious defect of PBE calculations in 
this system is that the energetics of the proton transfer 
is very poorly described as too much electronic charge 
is transferred with the proton. We show that a satis- 
factory description of the experimental data, including 
the absence of proton-transfer phases at low pressures, 
requires the inclusion of nuclear ZP motion and accurate 
descriptions, beyond those afforded by functionals such 
as PBE-GGA, of both exchange interactions and the van 
der Waals forces that arise from electron correlation. 



II. AB INITIO RANDOM STRUCTURE 
SEARCHING 

We have used the AIRSS method^ to identify low- 
enthalpy structures of AMH at pressures of up to about 
12 GPa. In the AIRSS approach randomly chosen struc- 
tures are relaxed to a minimum in the enthalpy at fixed 
pressure. In its simplest form AIRSS has almost no free 
parameters and is essentially unbiased, and it is there- 
fore the ideal basis upon which to impose constraints 
and biases towards the types of structure that one be- 
lieves are most favorable. Perhaps the simplest physical 
constraint that we have employed is to reject initial con- 
figurations in which atoms are closer than a defined mini- 
mum separation. One of the most useful constraints is to 
restrict the symmetries of the structures. The structures 
are chosen to obey the symmetries of a particular space 
group, although they are otherwise random, and the de- 
sired symmetry is maintained throughout the relaxation 



procedure. Another useful approach is to choose initial 
structures constructed from randomly placed "chemical 
units" , which in this case are equal numbers of NH3 and 
H2O molecules, or equal numbers of NH4 and OH units, 
or "hydrogen-bonded" NH 3 -H 2 AMH units. Each ini- 
tial unit cell was generated by choosing random unit cell 
translation vectors and scaling the volume to lie ran- 
domly within ±50% of some reasonable value. We then 
placed the required number of chemical units within the 
cell, applying symmetry constraints as required. Initial 
structures in which the overlap of molecules was signif- 
icant were rejected because they are likely to undergo 
unwanted chemical reactions during the relaxation pro- 
cedure. 

We first performed unconstrained searches with 2 for- 
mula units (f.u.), and we then searched with 4 f.u. and 
initial structures formed from NH3 and H2O molecules. 
Searches were then performed starting from random 
arrangements of either 2 or 4 preformed AMH units. 
Symmetry constrained searches were performed starting 
from random cells containing 2 randomly placed H2O 
molecules and 2 randomly placed NH3 molecules, and 
then applying the symmetry operations of space groups 
randomly chosen from those with 2 operations. Similar 
searches were performed with 4 symmetry operations and 
single units of H2O and NH3. To bias the procedure to- 
wards finding ionic structures, we started further searches 
with cells containing 4 f.u. and using building blocks of 
NH4 and OH units with n = 2 or n = 4 symmetry op- 
erations. 3 and 5 f.u. searches were performed without 
symmetry constraints. We also performed searches over 
larger unit cells containing 6, 6 and 8 f.u. in total, gener- 
ated with n — 6, n — 3 and n — 8 symmetry operations, 
respectively. A total of about 7,700 structures were re- 
laxed during the searches. 

The CASTEP plane wave code^ was used for all 
of the calculations on periodic crystals. Calculations 
were performed with the Perdew-Burke-Ernzerhof (PBE) 
generalized gradient approximation (GGA) exchange- 
correlation density functional^ the PBE functional 
with the Grimme semi-empirical dispersion correction 
(G06), 2? and the PBEO hybrid density functional^ which 
includes 25% exact (Hartree-Fock) exchange. We used 
ultrasoft pseudopotentials^i for the PBE and PBE+G06 
calculations and norm-conserving pseudopotentials gen- 
erated using the Opium software^ for the PBEO calcu- 
lations. For the searches we used a plane wave cut off 
energy of 340 eV and a Monkhorst-Pack ?6 Brillouin zone 
sampling grid of spacing 2n x 0.07 A -1 , and 2tt x 0.08 
A -1 for some of the searches with larger unit cells. All 
of the results reported in this paper were obtained by re- 
fining the structures obtained in the searches at a higher 
level of accuracy consisting of a plane wave cut off energy 
of 700 eV and a Brillouin zone sampling grid of spacing 
2-7T x 0.04 A -1 . The enthalpy difference between AMH 
I and the ionic PA/nmm phase of AMH reported here 
was changed by less than 0.0002 eV per AMH formula 
unit (f.u.) on doubling the cut off energy to 1400 eV, 



3 



0.8 



c 
3 



0.04 



^ 0.00 
-0.02 
-0.04, 
0.04 





AMH-I 
AMH-II (expt) 
AMH-II (calc) 


■ — ■ — P4/nmm (ionic) 
P2 1 2 1 2 1 (ionic) 






- a) 


1 ■ < 


1 ■ 1 
■ X' ' 


1,1,1,1, 







0.02 \ b) 



> 



0.00 
0.02 h 
0.04, 



2 

~T~ 



4 







0.04 



0.02- c ) 



§ 0.00 
-0.02 h 
-0.04, 



2 



4 

- r 



_l_ 



_i_ 



1 2 3 

Pressure (GPa) 



FIG. 1: Enthalpies per f.u. of AMH structures relative to 
AMH-II calculated with: a) the PBE functional, b) the PBE 
functional and G06 semi-empirical dispersion correction. The 
Gibbs free energy per f.u. of AMH structures relative to that 
of AMH-II is shown in c), calculated with the PBE functional 
including vibrational motion at 175 K. AMH-II (expt) de- 
notes the structure found in experiment while AMH-II (calc) 
denotes the structure found in the AIRSS study reported in 
Refi 



while the enthalpy change on doubling the number of 
k-points was even smaller. The required force tolerance 
for a successful geometry optimization in each search run 
was 0.05 eV/A, which was tightened to 0.005 eV/A for 
the final results reported in this paper. The stress on the 
unit cell was converged to better than 0.01 eV/A 3 . The 
norm-conserving pseudopotentials required a plane wave 
cut off energy of 1000 eV, for which the energy difference 
between AMH-I and PA/nmm was converged to better 
than 0.00002 eV per f.u. A less dense Brillouin zone sam- 
pling grid spacing of 2tt x 0.05 A- 1 was used for the PBE0 
calculations, which gave an energy convergence of better 
than 0.03 eV per f.u. 



III. RESULTS FROM STRUCTURE 
SEARCHING 

AIRSS was used successfully to determine the crystal 
structure of AMH-II in collaboration with experiment, 
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FIG. 2: Enthalpies per f.u. of PA/nmm relative to AMH- 
I for three methods; using the PBE functional, the PBE0 
hybrid functional, and PBE0 with both the G06 dispersion 
correction and ZP vibrational motion effects included. The 
PBE0 calculations were performed on the structures relaxed 
using the PBE functional. 
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FIG. 3: The PA/nmm ionic structure; light grey atoms are 
hydrogen, red atoms are oxygen, and blue atoms are nitrogen. 
The solid lines depict two unit-cells, the orientation of the 
crystallographic axes being indicated on the left. 



which provided initial constraints on the symmetry and 
dimensions of the unit ce\\X The AIRSS structure with 
112 atoms in the primitive unit cell was found to be al- 
most correct in a subsequent experiment, with the excep- 
tion that it produced one of the two possible H-bond or- 
dering schemes that is apparently not adopted by the real 
material. The relative enthalpies of the AMH-II struc- 
ture obtained from AIRSS and the experimental struc- 
ture, both fully relaxed within PBE, are shown in Fig. [1] 
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The PBE and PBE+G06 calculations indicate that the 
AMH-II structure obtained from AIRSS is roughly 0.01 
eV higher in enthalpy than the experimental structure. 
There is no reason why the experimental structure with 
the alternate H-bond ordering could not have been found 
if more searches had been performed. 

Fig.[lja) shows the variation of the enthalpy with pres- 
sure for a number of AMH phases calculated with the 
PBE functional. Our searches did not find any molecular 
structures more stable than AMH-I and II. This may sug- 
gest that, if any of the three unsolved AMH polymorphs 
are molecular, they are likely to have complex architec- 
tures with Z' > 8, where Z' is the number of molecules 
in the asymmetric unit cell. The ionic ammonium hy- 
droxide structure (space group P2\l\l\) was obtained in 
the DFT study of Fortes et al. 17 by compressing AMH-I 
which underwent an isosymmetric transition to the ionic 
(proton-transfer) form manifested by a discontinuity in 
the slope of the calculated energy- volume curve. This 
transition was also observed in the present work. The 
ionic P2i2i2i structure becomes more favourable than 
either AMH-I or AMH-II above about 6 GPa. Although 
the P2i2i2i ionic phase has a region of stability on this 
phase diagram relative to the known phases, there is no 
reason to believe that it is the most stable ionic phase 
of AMH, which motivates a systematic search, as de- 
scribed above. Searching with AIRSS revealed a struc- 
ture (shown in Fig. \S§ with space group PA/nmm to be 
the most stable in all 2 f.u. searches at 10 GPa. Sub- 
sequent searches with 4 f.u. at both 3 and 10 GPa also 
showed PA/nmm to be the most stable structure, re- 
gardless of the constraints imposed. Even searches with 
8 f.u. found PA/nmm to be the most stable. None of the 
searches performed with 3, 5, or 6 f.u. resulted in struc- 
tures with enthalpies as low as that of PA/nmm. In fact 
none of the space groups with 3 symmetry operations are 
subgroups of PA/nmm, and therefore it could not have 
been found in searches with 3 symmetry operations. 

The stability range of AMH-II between the AMH-I and 
ionic PA/nmm phases predicted by the PBE functional 
is very small (~0.01 GPa, see Fig. [1]). However, exper- 
iments have shown that the transition from AMH-I to 
AMH-II takes place at around 0.5 GPa, and then from 
AMH-II to AMH-IV (whose structure is presently un- 
known) at around 2.2 GPa at 170 KJ*2- An experimen- 
tal study has shown that warming AMH-IV at 6.5 GPa 
produces the body-centred-cubic (bec) phase VI at ^280 
Kji the structure of which has been reported to consist 
of orientationally and positionally disordered NH3 and 
H2O molecules.^ Our neutron powder diffraction study 
has shown that compression of AMH-V at room temper- 
ature does not lead to the formation of AMH- VI up to 
^9 GPa^i Interestingly, a similar disordered bec phase 
of ADH was reported by Fortes et al£2- and confirmed re- 
cently by Loveday et al^ It is likely that a solid solution 
could exist between the AMH and ADH compositions 
at high pressures, if the bec crystal structure is main- 
tained over a range of occupancies of the NH3 and H2O 



molecules. 

Clearly, construction of a complete computational 
phase diagram requires simulation of the AMH- VI struc- 
ture. However, the disordered nature of AMH- VI pre- 
cludes straightforward investigation using DFT but, as 
a first approximation, a so-called "shaking" search was 
performed, in which a 2x2x2 supercell comprising only 
the oxygen and nitrogen atoms was created. For each 
search the appropriate number of hydrogen atoms were 
distributed randomly over the supercell and the atomic 
positions were relaxed. The lowest enthalpy structure, 
which was obtained repeatedly, does not appear in Fig. 
[T]as it lies approximately 0.2 eV per f.u. above AMH-II 
and is quite far from thermodynamic stability. The large 
discrepancy between the computed stability and the re- 
producible experimental observation of AMH- VI is likely 
to be due to the small cell used in the calculations. 

PBE calculations predict the ionic PA/nmm phase to 
be thermodynamically stable across a broad region of the 
high pressure phase diagram. Whether PA/nmm corre- 
sponds to any of the polymorphs with as-yet undeter- 
mined structures (III, IV, or V), and its relationship (if 
any) to AMH- VI is not yet clear due to the lack of suit- 
able experimental data. The published neutron powder 
diffraction data for AMH-III, IV, and V are not of high 
quality?^ in particular the published powder pattern of 
AMH-IV consists only of very broad reflections. We have 
therefore carried out our own neutron powder diffraction 
studyj^i with the aim of obtaining high resolution data 
from these polymorphs. The results acquired so far con- 
firm that none of these polymorphs is likely to be the 
PA/nmm phase. There is, however, a discernible rela- 
tionship between the structures of the PA/nmm phase 
and AMH- VI, as described below. It is worth noting that 
apparently small structural changes, particularly those 
that break a crystal symmetry, can significantly affect 
a diffraction pattern. Therefore it is reasonable to ex- 
pect the diffraction patterns of PA/nmm and AMH- VI 
to differ, despite the relationship between the structures. 

With reference to Fig. SI the hydrogen-bonded layers 
in the PA/nmm ionic phase form a network with a square 
motif, which defines the tetragonal unit cell (marked in 
black). However, an oblique cell (dashed line) is also 
marked that closely approximates a cube having ammo- 
nium ions at the corners and a hydroxide ion near the 
center. A slight shrinkage of the tetragonal a- and b-axes 
(~4.4% relative to the value given in Table HIT) while keep- 
ing the length of the c-axis unchanged, forms a perfect 
cube. Furthermore, shifting the fractional z-coordinate 
of the oxygen atom from 0.6562 to 0.5 yields a heavy- 
atom structure of space-group Pm3m, a = 3.3850 A, 
and fractional atomic coordinates N = 0, 0, 0, and O 
= 0.5, 0.5, 0.5. Finally, mixing the occupancies of these 
two sites with ammonium and hydroxyl ions gives an 
ionic equivalent of the AMH- VI structure. In fact, the 
very small differences in computed Bragg intensities be- 
tween the ionic and molecular forms of AMH- VI lead us 
to conclude that AMH- VI may well be ionic rather than 
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FIG. 4: The ionic P4/nmm phase viewed along the c-axis. 
This structure comprises sheets of ions hydrogen-bonded to- 
gether (dashed rods) donated by the ammonium ions to the 
hydroxide. The hydrogen atom of the hydroxide ion points 
directly along the c-axis and does not appear to participate in 
a hydrogen bond. Hence, the layers are not H-bonded to one 
another and consequently the structure exhibits a large com- 
pressibility along the c-axis. The dashed line shows the basis 
of a cubic structure that may be derived from the tetragonal 
cell (marked by the solid black line) by a small distortion, as 
discussed in the text. 
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FIG. 5: Volume-pressure relationships for the AMH-I and 
AMH-II phases and the new ionic phase of AMH reported 
here, calculated with the PBE functional, including the ef- 
fects of ZP motion and temperature (175 K). The experi- 
mental data were obtained from deuterated samples at 180 K 
at the Institut Laue Langevin, in the gas phase with a high- 
resolution powder diffractometer, and from a Paris- Edinburgh 
cell study with a high-intensity powder diffractometer, for 
AMH-I and AMH-II, respectively. 



molecular, and that the PA/nmm phase might simply be 
an ordered variant. 



IV. ZERO-POINT MOTION AND FINITE 
TEMPERATURE EFFECTS 

The relatively small mass of hydrogen leads to large 
ZP motions in AMH, where 5 out of every 7 nuclei are 
protons. The ZP motion in AMH may lead to important 
differences in the relative stabilities of the phases, par- 
ticularly when comparing a dense ionic phase with a less 
dense molecular one. We have investigated the effects 
of ZP motion on the phase diagram of AMH within the 
quasi-harmonic approximation with the supercell method 
and finite atomic displacements. The quasi-harmonic ap- 
proximation normally gives a reasonable description of 
vibrational effects, including thermal expansion. 

We calculated the phonons of AMH-I and PA/nmm 
in 112 atom supercells, while a set of finite displacement 
phonon calculations were performed for AMH-II, which 
has a 112 atom primitive cell. 

Care was taken to ensure that the structures were very 
well relaxed prior to performing phonon calculations, en- 
suring that any stresses on the unit cell were less than 
0.01 eV/A 3 and that the forces were converged to within 
0.005 eV/A (0.003 eV/A for AMH-II). In addition, the 



fine grid on which the augmentation charge density for 
the ultrasoft pseudopotentials is representated was in- 
creased to 2.75 times the multiple of the wavefunction 
grid to obtain higher accuracy. 

The Gibbs free energy at 175 K is plotted against pres- 
sure for AMH-I, AMH-II and the P4/nmm ionic phase 
in Fig. [TJ as calculated with the PBE functional. The 
phonon pressure was evaluated from the derivative of 
the ZP energy (or Helmholtz free energy at finite tem- 
perature) with respect to volume. The total pressure is 
the sum of the static DFT and phonon pressures. The 
larger density of the ionic PA/nmm phase leads to higher 
phonon frequencies and hence destabilization relative to 
the molecular phases. The pressure obtained within PBE 
at which PA/nmm becomes the most stable at 175 K 
is increased by 0.55 GPa to 3.32 GPa. The inclusion 
of vibrational effects increases the pressure interval of 
stability of AMH-II to about 0.7 GPa, although this is 
still smaller than the experimental interval of about 1.7 
GPa.^ 

Experimentally, the transition from AMH-I to AMH-II 
at ^0.35 GPa results in a volume decrease of 4. 6%. 8 Both 
PBE, and PBE with ZP motion at 175 K, give only a 2% 
decrease in volume at this transition, whereas PBE with 
the G06 dispersion correction gives a decrease of 3.4%. 

Data from fits of the calculated pressure-volume data 
to the third-order Birch-Murnaghan equation of stated 
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9.7 

12.5 

12.2 

12.5 
8.9(4) 
7.33(3) 



5.0 
2.4 
2.4 
6.3 

4.2(3) 
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245.04 
241.50 
242.17 
218.10 
247.66 
248.00(2) 



8.1 
9.6 
9.2 
13.4 
7.2(3) 



5.5 
3.1 
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4.2 

5.3(2) 



971.60 
973.22 
983.25 
842.70 
947(2) 
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20.8 



13.8 
7.4 
5.4 

11.0 



106.51 
113.13 
117.95 
83.01 



Structure 



AMH-I (Z = 4) 



Method 



PBE 

PBE+ZP (0 K) 
PBE+ZP (175 K) 
PBE+G06 
Experiment (140 K) 1 
Experiment (180 K)£ 



AMH-II (Z = 16) 



PBE 

PBE+ZP (0 K) 
PBE+ZP (175 K) 
PBE+G06 
Experiment (180 K) 



PA/nmm (Z = 2) 



PBE 

PBE+ZP (0 K) 
PBE+ZP (175 K) 
PBE+G06 



TABLE I: The bulk modulus (K), the first pressure derivative of the bulk modulus (K'), and the equilibrium volume (Vo). 
Fitting ranges: AMH-I 0-2.5 GPa, AMH-II 0.4-2.5 GPa. PA/nrnrn PBE 3-6 GPa, PBE+G06 0-2.5 GPa. 



Functional 


Lattice parameters 

(A,°) 




Atom 


Site 


Wyckoff 
symbol 


Fractional atomic coordinates 


PBE 


a=5.006 6=5.006 c= 


=3.385 


N 


-42 


2a 


0.25 


0.7500 


0.0000 




q=90 /3=90 7= 


=90 


O 


4mm 


2c 


0.25 


0.2500 


0.6562 








HI 


4mm 


2c 


0.75 


0.7500 


0.6308 








H2 


.m. 


8i 


0.25 


0.5727 


0.8262 



TABLE II: Structure of the P^jnmm ionic phase at 3 GPa with the PBE functional (primitive cell, Z — 2). 



with and without the ZP motion and temperature con- 
tributions, are shown in Table QJ The uncorrected PBE 
results compare most favorably with the experimental 
equilibrium volume, bulk modulus and the first pres- 
sure derivative of the bulk modulus for both AMH-I and 
AMH-II. It is interesting to note that there is one excep- 
tion to the expected increase in volume from including 
ZP motion, AMH-I is seen to shrink slightly when ZP 
motion is taken into account. All three structures are 
found to increase in volume when thermal effects are in- 
cluded at 175 K, accompanied by a small reduction in the 
bulk moduli relative to the values obtained on including 
ZP motion at K. 




2 3 4 

Pressure (GPa) 



FIG. 6: Variation of the lattice parameters with pressure for 
the P4/nmm ionic phase, calculated with the PBE functional, 
with and without the G06 dispersion correction. 



V. DISPERSION CORRECTION 

The transition pressure between the AMH-I and II 
molecular phases obtained with PBE of 2.8 GPa is sig- 
nificantly larger than the experimental value of 0.5 GPa. 
Here we explore the effects of including dispersion forces 



which are not described by density functionals such as 
PBE. For this purpose we have recalculated the phase di- 
agram using the PBE functional with the Grimme semi- 
empirical dispersion correction (G06)^ see Fig. [T] b). 
The significant overestimate of the transition pressure 
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Method 


AMH I -> AMH II 


Transition Pressure (GPa) 
AMH II PA/nmm 


AMH I -> PA/nmm 


PBE 


2.7 


2.8 


2.8 f 


PBE+G06 


0.9 f 


0.4 f 


0.5 


PBE+ZP (175 K) 


2.7 


3.3 


3.3 f 


PBEO 


n/a 


n/a 


10.8 


PBE0+ZP+G06 (0 K) 


n/a 


n/a 


8.8 


Experiment (180 K)I 


0.5 


n/a 


n/a 



TABLE III: Summary of the transition pressures between the AMH phases studied here, with each of the methods utilised. 
Transitions which would occur at a pressure where neither phase is thermodynamically stable are denoted by Note that 
PBEO calculations were not performed for the AMH II phase and are therefore not included in the table (n/a). 



obtained with PBE is substantially improved by includ- 
ing the G06 correction. The transition pressure between 
the molecular AMH-I and AMH-II phases is reduced to 1 
GPa using the PBE+G06 functional, which is still some- 
what larger than the experimental value of 0.5 GPa. 
However, the G06 correction significantly favours the 
denser PA/nmm ionic phase, which now becomes stable 
at about 0.5 GPa. Each phase undergoes a substantial 
volume contraction when the G06 correction is included, 
which is likely due to an overestimation of dispersion ef- 
fects. 

The large effect of the dispersion correction on the 
ionic structure may seem paradoxical as it amounts to 
a relatively small fraction of the binding energy of the 
NH|" ■ • -OH - ionic complex (see Fig. [8]). The rapid re- 
duction in the volume of PA/nmm with applied pressure 
apparent in Fig. [5] is almost entirely associated with a 
contraction along the c lattice parameter. As shown in 
Fig. HI the hydrogen bonds form a square net within the 
a-b-planes of PA/nmm, but there is no apparent hydro- 
gen bonding between the layers (see Fig. [3]), and therefore 
it is very soft in the c direction. 



VI. BEYOND THE PBE FUNCTIONAL 
A. Calculations for molecular and ionic fragments 

Despite giving a good description of the bulk structural 
properties of AMH-I and II, the PBE functional leads to 
substantially incorrect transition pressures. PBE overes- 
timates the pressure of the AMH-I/II transition, and the 
PA/nmm ionic phase is predicted to be sufficiently stable 
to virtually eliminate the region of stability of AMH-II. 
Even accounting for ZP motion and thermal effects it is 
hard to reconcile the apparent stability of the PA/nmm 
ionic phase with the existing room temperature experi- 
ments in which the ionic AMH- VI phase was not found 
at pressures as high as 9 GPa^I The structures of AMH- 
III, IV and V have not yet been solved, and one or more 
of these phases may be ionic. 

To obtain insight into this problem we have selected 



fragments of the AMH-I and PA/nmm crystals (at pres- 
sures of 1 and 4 GPa, respectively) for a more detailed 
analysis. These fragments were chosen to be represen- 
tative of the interactions in the crystal and consisted of 
pairs of molecules/ions in close proximity: NH 3 and H 2 
for the AHM-I crystal (Fig. ED and NHj and OH' ions 
for the PA/nmm crystal (Fig. [5]). In Figs. [7] and [5] we 
plot interaction energies for these pairs calculated using 
the supermolecular approach, that is, Ei n % = E(AB) — 
E(A) - E{B), where E(A) and E{B) are the energies of 
the monomers and E(AB) is the energy of the complex. 
We have used the CCSD(T), MP2, PBE, PBE+G06 and 
PBEO methods, though the MP2 results are not shown as 
they are essentially identical to the CCSD(T) values. We 
have calculated counterpoise-corrected interaction ener- 
gies using the Boys and Bcrnardi 32 scheme with a Sadlej- 
pVTZ basis^ augmented with a small set of "bond cen- 
tered functions"— which help to saturate the dispersion 
energy^ Selected calculations with the larger aug-cc- 
pVTZ basis set suggest that the CCSD(T) interaction 
energies are converged to 5% at the equilibrium geome- 
try and better at larger separations. These calculations 
were performed using the DALTON 2.0 program^ For 
deeper insight into the nature of the interaction energies 
we have additionally used the CamCASP program^! to 
perform symmetry adapted perturbation theory (SAPT) 
DFT— calculations to decompose the interaction en- 
ergies into physical components. 

Fig.[7]shows that PBE overbinds the NH 3 - • -H 2 com- 
plex, although the equilibrium bond length agrees reason- 
ably well with the CCSD(T) value. The overestimation of 
the bulk modulus (K) of AMH-I reported in Table U may 
be a reflection of the overbinding of PBE and the conse- 
quent overestimate of the curvature of the potential well. 
In contrast, PBEO underbinds the complex, but results in 
an equilibrium separation in near perfect agreement with 
the CCSD(T) calculations. The SAPT(DFT) energy de- 
composition shows that asymptotically the NH3- • -E^O 
interaction is dominated by the dipole-dipole electro- 
static energies with the dispersion and polarization en- 
ergies being negligible. Consequently it should not be a 
surprise that both the correlated and density-functional 
methods agree in this region. 
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FIG. 7: Interaction energy against separation of NH3 and 
H2O for a configuration taken from the AMH-I molecular 
phase. 
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FIG. 8: Interaction energy against separation of NHj and 
OH~ for a configuration taken from the PA/nmm ionic phase. 



The ionic system behaves very differently from the 
molecular one. Note that the interaction energies for the 
ionic system shown in Fig. [8] are an order of magnitude 
larger than in the molecular system. The SAPT(DFT) 
results show that the energetics are almost completely 
dominated by the charge-charge electrostatic interaction 
between the ions, with the dispersion energy contribu- 
tion being insignificant even at the relatively short N- ■ O 
separation of 2.5 A. This is exactly the kind of system 
for which local and semi-local density functionals are ex- 



pected to achieve high accuracy. Indeed, the PBE and 
PBEO energies are close to the CCSD(T) result at the 
equilibrium geometry and on the repulsive wall. In con- 
trast to the molecular system, however, PBE significantly 
overbinds the ionic complex at large separations. As 
can be seen in Fig. [SJ the overbinding compared with 
CCSDfT) is 20% at 5 A but it grows to more than 50% 
at 6.5 A. The interaction energies are still relatively large 
at these separations, and consequently these errors have 
important effects within the crystal. This is consistent 
with our observation that the PA/nmm ionic phase is 
over-stabilized with respect to the molecular phases. 

A partial charge analysis shows that the substantial 
difference between the PBE and CCSD(T) interaction 
energies at large separations is a consequence of exces- 
sive charge transfer. Within PBE the magnitude of the 
charges on the monomers grows as they are pulled apart, 
becoming as large as ±1.2e at a N- ■ O separation of 6.5 
A. In contrast, the CCSD(T) potential energy curve can 
be fitted very well with a —1/R form corresponding to 
partial charges of ±le. The PBE charge transfer error is 
analogous to the derealization or static correlation errors 
exhibited by local and semi-local density functionals. 41 
This error can be at least partially corrected by intro- 
ducing some fraction of non-local exchange. Using PBEO 
gives a significant improvement upon the PBE energies, 
with the overbinding compared with CCSD(T) at 5 A 
and 6.5 A reduced to 6% and 28%, respectively, see Fig. 

E 



B. PBEO calculations for crystalline AMH 

As the PBEO functional appears to offer a partial solu- 
tion to the charge transfer errors which arise with PBE, 
we performed calculations for the molecular AMH-I and 
ionic PA/nmm phases using the PBEO hybrid functional. 
Because of the slow convergence of the exchange terms 
with distance, the PBEO calculations for the crystal were 
as much as two orders of magnitude more computation- 
ally expensive than the corresponding PBE calculations. 
For this reason geometry optimizations were not possible 
with PBEO and instead we performed single-point energy 
calculations using the relaxed PBE structures. We at- 
tempted to obtain equation of state parameters from the 
PBEO results for AMH-I and PA/nmm. The calculated 
energies were fitted to the third-order Birch-Murnaghan 
equation of state and the pressure was obtained by differ- 
entiation. We obtained equilibrium volumes of Vq = 220 
A 3 for AMH-I and V = 172 A 3 for PA/ nmm. which 
are similar to those obtained from the PBE+G06 calcu- 
lations. We were not able to obtain reliable estimates of 
K or K' from our PBEO data as the values were sensitive 
to the data and fitting procedure. The use of the PBE 
structures for the PBEO calculations is an approximation. 
However, we believe this approach to be reasonably ro- 
bust because enthalpies obtained using partially relaxed 
PBE structures were almost identical to those from the 
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FIG. 9: Variation of the band gaps of the phases with pres- 
sure. AMH-I (and AMH-II using PBE only) have direct band 
gaps, while both the PBE and PBEO calculations predict that 
the PA/ nmm phase has an indirect band gap. The PBEO cal- 
culations were performed using the PBE structures. 



fully relaxed structures. 

Mulliken charge analysis of the PA/ nmm crystal sug- 
gests that the electronic charge transfer between the 
"OH"" and "NIT|" ions has a magnitude of about 0.70e 
when using the PBE functional and 0.55e with PBEO. 
The improved description of the ionic phase provided by 
PBEO destabilizes the PA/ nmm ionic phase relative to 
AMH-I, and consequently PA/ nmm enters above 10.8 
GPa, as can be seen in Fig. H AMH-II, AMH-IV and 
AMH-VI are not present on this phase diagram, but they 
are expected to have regions of stability between AMH-I 
and PA/nmm, pushing the entry of the new ionic phase 
to even higher pressures. The inclusion of vibrational ef- 
fects further destabilizes PA/nmm and it only becomes 
stable at 11.3 GPa, although conversely the inclusion of 
the G06 dispersion correction reduces the transition pres- 
sure to 8.8 GPa, as shown in Fig. [5] 



VII. BAND GAPS OF THE PHASES 



PBE functional. A new ionic structure of space group 
PA/nmm was found to be stable above 2.8 GPa. Subse- 
quent investigations into the effects of temperature, ZP 
motion and dispersion forces found that the latter two 
play a substantial role in determining the relative stabil- 
ities of the phases. The inclusion of ZP motion desta- 
bilizes the dense PA/nmm ionic phase; conversely the 
dispersion correction leads to a significant underestima- 
tion of the transition pressure from molecular AMH to 
ionic PA/nmm. The dispersion forces are attractive and 
therefore tend to favour denser structures. 

The relationship described earlier between the struc- 
tures of the PA/nmm phase and AMH-VI may indi- 
cate that we have discovered an ordered ionic variant 
of AMH-VI. Alternatively, the PA/nmm phase may be 
related to one of the three experimentally observed AMH 
phases (III, IV, V) whose structures are unknown. As the 
PA / nmm phase was the only new structure found which 
is predicted to be thermodynamically stable within some 
pressure range, it is likely that if it is one of the three 
unknown AMH structures, then the other two will have 
unit cells with Z > 8. 

Accurate CCSD(T) calculations on representative 
complexes from the ionic and molecular AMH phases re- 
vealed that the PBE functional substantially over/binds 
the ionic phase. We have presented evidence that this 
overbinding arises from an overestimate of the electronic 
charge transfer accompanying the proton transfer, which 
can be partially remedied by using the hybrid PBEO 
functional. Using the PBEO functional leads to an in- 
crease in the enthalpy of the ionic PA/nmm phase by 
about 0.6 eV per f.u. relative to AMH-I. The transition 
pressure from AMH-I to PA/ nmm phase is substantially 
increased, which eliminates the inconsistency with ex- 
periment. This failure of the PBE functional, for an 
electrostatic-bound system for which GGA-type density 
functionals are typically assumed to be accurate, is likely 
to have implications for other systems. Further experi- 
mental work is necessary to explore the phase diagram 
of AMH at 9 GPa and above to confirm our assignment 
of PA/ nmm as a stable ionic phase of AMH. 



The minimum band gaps obtained with the PBE and 
PBEO functionals are shown in Fig. [^calculated using the 
LinDOS code^S While the molecular structures have di- 
rect band gaps, both PBE and PBEO calculations predict 
PA/ nmm to have an indirect band gap. The band gaps 
provided by the PBEO functional, which includes explicit 
exchange interactions, are substantially larger than those 
predicted by PBE. 

VIII. CONCLUSIONS 

Wc have performed extensive ah initio searches for new 
phases of AMH at pressures up to 12 GPa using the 
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